#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fit.h>

#define RAW(v) ((void*)&(v))

void test_gsl_linear_regression(void)
{
    int i = 0;

    //sample
    size_t norris_n = 11;
    size_t xstride = 1;
    size_t ystride = 1;

    double x[] = { 71,  73,  64,  65,  61,  70,  65,  72,  63,  67,  64};
    double y[] = {160, 183, 154, 168, 159, 180, 145, 210, 132, 168, 141};

    //result
    double c0 = 0.0;
    double c1 = 0.0;
    double cov00 = 0.0;
    double cov01 = 0.0;
    double cov11 = 0.0;
    double sumsq = 0.0;

    //test infinite.
    //for (i = 0;i < norris_n;i++) {
    //    x[i] = 10;
    //}

    //calculate
    gsl_fit_linear (x, xstride, y, ystride, norris_n, 
            &c0, &c1, &cov00, &cov01, &cov11, &sumsq);

    if (gsl_finite(c0)) {
        printf("f(x) = %lf + %lf x\n", c0, c1);
        printf("cov00 = %lf\n", cov00);
        printf("cov01 = %lf\n", cov01);
        printf("cov11 = %lf\n", cov11);
        return;
    }

    //exchange and try again
    gsl_fit_linear (y, ystride, x, xstride, norris_n, 
            &c0, &c1, &cov00, &cov01, &cov11, &sumsq);

    if (gsl_finite(c0)) {
        printf("f(y) = %lf + %lf y\n", c0, c1);
        printf("cov00 = %lf\n", cov00);
        printf("cov01 = %lf\n", cov01);
        printf("cov11 = %lf\n", cov11);
        return;
    }

    //error
    printf("Error!\n");
    return;
}

void test_gsl_lu_decomp(void)
{
    //Ax=b
    gsl_matrix* A = gsl_matrix_alloc(2,2);
    gsl_matrix* a = gsl_matrix_alloc(2,2);
    gsl_vector* b = gsl_vector_alloc(2);
    //result
	gsl_permutation* p = gsl_permutation_alloc(2);
	gsl_vector *sx = gsl_vector_alloc(2);
    int sig = 0;

    int i = 0;
    int j = 0;

    //Clear data
    gsl_matrix_set_zero(A);
    gsl_vector_set_zero(b);

    //2x + y  = 1
    //2x + 6y = -9
    //
    //(x, y) = (3/2, -2)
    gsl_matrix_set(A, 0, 0, 2.0);
    gsl_matrix_set(A, 0, 1, 1.0);

    gsl_matrix_set(A, 1, 0, 2.0);
    gsl_matrix_set(A, 1, 1, 6.0);

    gsl_vector_set(b, 0, 1.0);
    gsl_vector_set(b, 1, -9.0);

    //Decompress PA = LU
    //so: PAx = (LU)(x) = Pb
    //  Ly = Pb
    //  Ux = y
    //this will modify A! so duplicate A as a.
    //(A, p) -> (LU, p, sig)
    gsl_matrix_memcpy(a, A);
	gsl_linalg_LU_decomp(a, p, &sig);
	gsl_linalg_LU_solve(a, p, b, sx);

    printf("Ax = b\n");
    for (i = 0;i < 2;i++) {
        printf("|");
        for (j = 0;j < 2;j++) {
            printf("%6.2lf ", gsl_matrix_get(A, i, j));
        }
        printf("|");
        printf(" %6.2lf ", gsl_vector_get(sx, i));
        printf("|");
        printf(" %6.2lf", gsl_vector_get(b, i));
        printf("|");
        printf("\n");
    }

    gsl_matrix_free(A);
    gsl_matrix_free(a);
    gsl_vector_free(b);
    gsl_permutation_free(p);
    gsl_vector_free(sx);
}

int main(int argc, char* argv[])
{
    test_gsl_linear_regression();
    test_gsl_lu_decomp();
    return 0;
}

